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ABSTRACT 

In  the  present  investigation  both  2D  and  3D  linear  elastic  analyses  have  been 
undertaken  to  assess  the  accuracy  and  computational  resources  associated  with  the 
use  of  h-element  and  p-element  formulations.  All  analyses  have  been  undertaken  for 
a  local  stress-concentrating  feature,  which  is  typical  of  aircraft  structures,  namely  a 
circular  hole  in  a  remotely  loaded  plate.  Here  highly  accurate  results  are  obtained  for 
both  element  formulations,  and  it  is  found  that  die  p-elements  offer  large  savings  in 
analysis  times.  Subsequently  a  relationship  is  developed  and  validated  for  both  2D 
and  3D  meshes  to  determine  an  equivalent  non-uniform  h-element  mesh  density, 
which  will  yield  the  same  accuracy  as  a  fully  converged  p-element  mesh.  This 
provides  a  useful  transferable  tool  for  designing  cost-effective  h-element  meshes, 
obviating  the  need  for  multiple  mesh  refinement  iterations.  Finally,  it  is 
demonstrated  numerically  that  p-elements  of  order  five  and  above,  can  be  used  to 
predict  accurate  through-thickness  peak  stresses,  in  a  typical  plate  modelled  with 
only  one  element  through  its  full  thickness. 
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Executive  Summary 


Finite  element  methods  (FEM)  are  currently  used  extensively  throughout  the  Airframes 
and  Engines  Division  (AED)  as  part  of  its  role  in  providing  through  life  support  of 
aircraft  in  service  with  the  Royal  Australian  Air  Force  (RAAF).  In  particular,  stress 
predictions  obtained  from  FEM  are  used  in  designing  and  validating  life  extension 
repairs,  assessing  the  structural  integrity  of  aircraft,  and  in  the  life  assessment  of  aging 
aircraft.  Typically,  the  role  of  the  FE  analysis  is  to  determine  accurate  numerical 
predictions  of  peak  stresses  at  local  geometric  features  such  as  notches,  in  large  and 
complex  components,  such  as  the  wing  pivot  fitting  of  the  F-lll  aircraft  in  service  with 
the  RAAF.  Historically,  the  commercial  finite  element  codes  used  in  AED  to  undertake 
such  stress  analyses  have  used  h-element  formulations.  In  recent  years  p-element 
formulations  have  become  available  in  some  commercial  codes,  which  potentially  offer 
significant  advantages  in  providing  more  efficient  and  cost  effective  solutions. 
However,  there  appears  to  be  limited  published  results  available  for  2D  and  3D  notch 
stress  evaluation  using  p-elements. 

Hence,  in  the  present  investigation  both  2D  and  3D  linear  elastic  analyses  have  been 
undertaken  to  assess  the  accuracy  and  computational  resources  associated  with  the  use 
of  h-element  and  p-element  formulations.  All  analyses  have  been  undertaken  for  a 
local  stress-concentrating  feature,  which  is  typical  of  aircraft  structures,  namely  a 
circular  hole  in  a  remotely  loaded  plate.  Here  highly  accurate  results  are  obtained  for 
both  element  formulations,  and  it  is  found  that  the  p-elements  offer  large  savings  in 
analysis  times.  Subsequently  a  relationship  is  developed  and  validated  for  both  2D  and 
3D  meshes  to  determine  an  equivalent  non-uniform  h-element  mesh  density,  which 
will  yield  the  same  accuracy  as  a  fully  converged  p-element  mesh.  This  provides  a 
useful  transferable  tool  for  designing  cost-effective  h-element  meshes,  obviating  the 
need  for  multiple  mesh  refinement  iterations.  Finally,  it  is  demonstrated  numerically 
that  p-elements  of  order  five  and  above,  can  be  used  to  predict  accurate  through¬ 
thickness  peak  stresses,  in  a  typical  plate  modelled  with  only  one  element  through  its 
full  thickness. 

The  results  obtained  in  the  present  work  provide  valuable  transferable  outputs,  which 
can  be  used  to  undertake  cost  effective  analyses  for  both  h-element,  and  p-element 
technologies,  relating  to  the  structural  assessment  of  RAAF  aircraft.  It  is  noteworthy 
that  very  recently  these  results  have  been  used  for  efficient  FE  analysis  of  large-scale 
structural  models,  as  part  of  the  Australian  F-lll  Aircraft  Structural  Integrity  Sole 
Operator  Program. 
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1.  Introduction 


Finite  element  methods  (FEM)  are  currently  used  extensively  throughout  the  Airframes 
and  Engines  Division  (AED)  as  part  of  its  role  in  providing  through  life  support  of 
aircraft  in  service  with  the  Royal  Australian  Air  Force  (RAAF).  In  particular,  stress 
predictions  obtained  from  FEM  are  used  in  designing  and  validating  life  extension 
repairs,  for  assessing  the  structural  integrity  of  aircraft,  and  in  the  life  assessment  of 
ageing  aircraft.  Typically,  the  role  of  the  FE  analysis  is  to  determine  accurate  numerical 
predictions  of  peak  stresses  at  local  geometric  features  such  as  notches  in  large  and 
complex  components,  such  as  the  wing  pivot  fitting  of  the  F-lll  aircraft  in  service  with 
die  RAAF.  Historically,  the  commercial  finite  element  codes  used  in  AED  to  undertake 
most  of  the  stress  analyses  have  used  h-formulation  elements1.  In  recent  years, 
p-formulation  elements1  have  become  available  in  some  codes,  which  potentially  offer 
significant  advantages  in  providing  more  efficient  and  cost  effective  solutions. 

In  the  FEM  of  computational  stress  analysis,  a  component  is  modelled  as  an 
assemblage  of  discrete  (or  finite)  elements.  If  conventional  h-element  FE  technology  is 
adopted  [1],  the  analysis  error  in  predicted  peak  stresses  generally  decreases  as  the  size 
of  elements  decreases.  This  is  because  typically  each  element  can  accurately  model  a 
linear  variation  in  stress,  however  throughout  the  structure  the  stress  distribution  is 
significantly  non-linear.  Hence  to  obtain  accurate  stress  results,  a  modeller  would 
ideally  perform  several  analyses  of  the  same  geometric  model  using  an  increasing 
number  of  smaller  elements,  until  convergence  was  achieved  (i.e.  stopping  when  the 
peak  stress  difference  between  successive  analyses  falls  to  an  acceptable  level).  One 
difficulty  with  this  approach  is  that  a  new  geometric  element  mesh  must  be  created  for 
each  successive  analysis.  Remeshing  can  be  a  highly  time  consuming  task  for  the 
modeller,  even  if  state-of-the-art  pre-processing  tools  are  used  (particularly  for  models 
which  have  cracks  or  are  three-dimensional).  Also,  progressively  greater 
computational  resources  are  then  needed  in  the  solution  phase. 

Unlike  h-elements,  p-elements  are  defined  by  variable  order  polynomial  displacement 
functions,  and  hence  can  potentially  model  higher  order  stress  variations  in  each 
element  accurately.  Here  to  achieve  convergence,  instead  of  reducing  the  element  size, 
only  the  polynomial  order  of  the  elements  needs  to  be  increased.  However,  while  p- 
elements  have  been  reported  in  the  literature  for  several  years  [2-5],  it  appears  that 
comparisons  of  the  relative  performance  of  p-elements  and  h-elements  have  only  been 
made  for  some  limited  two-dimensional  (2D)  problems  [6].  Also,  there  appear  to  be 
limited  publications  relating  to  the  use  of  p-elements  for  more  complex  applications, 
and  only  in  recent  years  have  p-elements  become  available  to  some  extent  in 
commercial  FE  codes.  For  instance,  p-elements  are  available  in  MSC.NASTRAN  [7]  for 
linear  elastic  analysis,  but  not  for  plasticity  analysis,  and  ABAQUS  [8]  and  PAFEC  [9] 
do  not  have  p-elements  at  all.  Hence,  before  the  p-element  approach  can  be  used  with 

1  For  brevity,  h-formulation  elements  will  be  referred  to  as  h-elements  and  p-formulation 
elements  as  p-elements  in  the  rest  of  the  report. 


1 


DSTO-TR-1151 


confidence  in  AED,  it  is  considered  beneficial  to  undertake  some  appropriate 
benchmark  studies  comparing  the  performance  of  the  two  element  types  for  the 
calculation  of  typical  notch  stresses,  in  two  and  three  dimensions.  Also  of  further 
interest  is  the  hypothesis  that  even  if  a  FE  code  does  not  support  p-element  technology, 
an  alternate  p-element  analysis  in  another  code  may  be  useful  in  guiding  the  design  of 
an  equivalent  h-element  mesh,  without  the  need  for  time  consuming  multiple  mesh 
refinement  iterations. 

Hence  in  the  present  investigation,  both  2D  and  3D  analyses  were  undertaken  to  assess 
the  accuracy  and  computational  resources  needed  for  various  uniformly  distributed 
h-element  and  p-element  mesh  densities,  and  p-element  orders.  All  analyses  are 
undertaken  for  a  local  stress-concentrating  feature,  which  is  typical  of  aircraft 
structures,  namely  a  circular  hole  in  a  remotely  loaded  plate.  From  these  results,  an 
empirical  relationship  is  then  developed  for  both  2D  and  3D  meshes  to  determine,  in 
one  iteration,  the  non-uniform  h-element  mesh  density  which  will  yield  approximately 
the  same  accuracy  as  a  given  non-uniform  order  p-element  mesh.  Through  a  numerical 
example,  the  veracity  of  this  approach  is  then  investigated  as  a  potential  transferable 
tool  in  designing  cost-effective  h-element  meshes.  Finally,  an  investigation  is 
undertaken  to  determine  the  p-element  modelling  requirements  for  the  accurate  and 
efficient  prediction  of  peak  (and  non-uniform)  through-thickness  stresses  in  plates. 


2.  Geometry  and  Finite  Element  Approach 


Circular  holes  in  loaded  plates  are  a  common  stress-concentrating  feature  of  aircraft 
structures,  and  hence  this  is  the  case  that  has  been  chosen  for  the  present  numerical 
benchmark  analyses.  The  geometry  under  consideration  is  shown  in  Figure  1.  Here  the 
flat  plate  contains  a  centrally  located  hole  of  radius  r  and  is  loaded  by  a  uniaxial  remote 
stress  Si ,  on  the  faces  at  x  =  ±  10.177r,  in  the  x  direction.  The  plate  width  is  8.496r  while 
the  thickness  is  1.77r.  This  specific  plate  geometry  has  been  chosen  to  allow  convenient 
comparison  with  prior  published  3D  FE  results  [10].  In  the  present  work  one  plane  of 
symmetry  was  invoked  for  all  2D  analyses,  while  typically  three  planes  of  symmetry 
were  utilised  for  the  3D  analyses.  For  all  analyses  the  hole  radius  was  taken  as 
r  =  11.3  mm  and  remote  stress  was  Si  =  30  MPa.  Also,  the  numerical  values  for  the 
material  constants  were  chosen  to  be  representative  of  an  aluminium  alloy,  with 
Young's  Modulus  taken  as  73049  MPa,  and  Poisson's  ratio  taken  as  0.33. 

All  analyses  were  undertaken  using  a  Hewlett-Packard  K  Series  9000  computer  located 
in  AED.  The  finite  element  codes  PAFEC  (Version  8.5)  and  MSC.NASTRAN  (Version 
69)  were  used  for  the  h-element  runs,  and  the  p-element  runs  were  performed  in 
MSC.NASTRAN  (Versions  69  and  70.5).  Pre-processing  and  post-processing  for  the 
PAFEC  and  MSC.NASTRAN  analyses  were  undertaken  using  PIGS  and 
MSC.PATRAN  (Versions  7.0A  and  7.5)  respectively. 
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For  the  p-element  analyses  it  was  important  to  remove  all  possible  sources  of 
singularity  for  this  study,  because  very  accurate  results  were  required.  Hence  typically 
the  following  techniques  were  adopted  in  the  modelling  process:  (i)  loads  were  applied 
as  pressures  on  element  edges  (2D  cases)  or  faces  (3D  cases)  rather  than  as  forces  on 
nodes,  and  similarly  (ii)  constraints  were  applied  to  element  edges  or  faces  rather  than 
to  nodes.  Also,  a H  element  edges  are  considered  to  be  straight  unless  specified 
otherwise.  It  is  noted  that  when  using  MSC.NASTRAN,  to  obtain  highly  accurate 
results,  the  shape  of  curved  external  edges  must  be  specified  precisely.  Hence  for  the 
analyses  presented  in  Sections  3  and  4  the  edges  were  defined  by  equations 
representing  a  circular  arc.  This  method  required  manual  insertion  of  equations  and 
element  edge  definitions  into  the  MSC.NASTRAN  input  deck.  However,  based  on  the 
work  in  Sections  3  and  4,  it  was  concluded  that  this  approach  would  be  too  time 
consuming  for  the  modelling  of  complex  structures.  Hence,  the  final  set  of  3D  analyses 
(see  Section  5)  was  performed  with  32-noded  brick  elements,  which  have  two  midside 
nodes  on  each  element  edge  [7],  These  midside  nodes  were  used  to  implement  the  two- 
point  method  for  defining  the  curvature  of  curved  element  edges  and  faces  in 
MSC.NASTRAN. 

An  error  measure  for  predicted  Kt  was  used  to  assess  the  accuracy  of  the  various 
modelling  techniques  and  mesh  densities.  In  this  study  the  error  was  defined  as: 


SKt  ~  (Kt(FEA)  "  Kt(exact))/  Kt(exact)  (1) 

where  Kt(FHA)  is  the  maximum  stress  concentration  factor  from  the  analysis  under 
consideration,  and  Kt(exact)  is  the  exact  solution  for  maximum  stress  concentration  factor. 
The  present  geometry  is  a  finite  width  plate,  for  which  no  closed  form  theoretical 
solution  appears  to  be  available,  although  an  approximate  solution  obtained  by 
Howland  gives  a  value  of  Kt=  3.21  [11].  Hence  in  the  present  work,  the  'exact'  solution 
is  taken  to  be  the  result  from  the  most  converged  analysis,  which,  for  both  the  2D  and 
3D  benchmarks,  was  obtained  from  a  high  order  p-element  analysis. 


3.  Convergence  of  2D  P-Element  and  H-Element 

Benchmarks 


3.1  Convergence  Analyses 

In  this  section  benchmark  results  are  given  for  several  variations  of  2D  p-element  and 
h-element  meshes.  In  these  plane-stress  analyses  the  plate  thickness  was  arbitrarily 
selected  as  10  mm.  MSC.NASTRAN  was  used  for  the  p-element  runs  and  PAFEC  for 
the  h-element  analyses,  with  the  aim  of  comparing  the  accuracy  and  efficiency  of  the 
two  element  formulations.  In  all  analyses  the  element  mesh  density  was  approximately 
uniform.  The  PAFEC  h-element  analyses  used  8-noded  iso-parametric  quadrilateral 
elements,  whereas  4-noded  quadrilateral  p-elements  (which  can  have  curved  edges) 
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were  used  in  the  MSC.NASTRAN  analyses.  For  the  PAFEC  h-element  analyses,  the 
mesh  density  was  varied  from  one  division  per  half-width  of  the  plate  to  32  divisions 
per  half-width,  as  shown  in  Figure  2.  Hence,  the  number  of  h-elements  varied  from  10 
to  10240.  The  base  model  for  the  MSC.NASTRAN  p-element  analyses  contained  40 
elements  as  shown  in  Figure  2(b),  i.e.  the  same  mesh  as  the  PAFEC  h-element  analysis 
with  two  divisions  per  half-width.  Seven  p-element  analyses  were  conducted  with  the 
polynomial  order  of  the  elements  taking  values  from  one  to  seven  [12],  with  all 
elements  having  the  same  polynomial  order  in  an  individual  analysis. 

3.2  Convergence  Results 

For  all  analyses,  results  for  the  following  were  obtained:  (i)  number  of  degrees  of 
freedom2,  (ii)  total  analysis  time,  (iii)  maximum  stress  concentration  factor,  and 
(iv)  percent  error  in  the  maximum  stress  concentration  factor.  These  results  are  listed  in 
Tables  1  and  2  for  the  h-element  and  p-element  analyses  respectively.  Table  2  also  gives 
results  obtained  from  the  MSC.NASTRAN  error  estimator  at  the  grid  point  (i.e.  node) 
location  with  the  maximum  stress  concentration  factor.  The  MSC.NASTRAN  grid  point 
error  estimator  [7]  is  based  on  the  differences  in  stress  computed  from  all  elements  that 
use  the  designated  grid  point  (i.e.  grid  point  stress  discontinuity).  Because  a  converged 
or  'exact'  value  of  maximum  stress  concentration  factor  is  available,  it  is  possible  to 
assess  the  accuracy  of  the  MSC.NASTRAN  error  estimator  for  p-elements.  As  expected, 
it  was  found  that  the  maximum  stress  concentration  always  occurred  at  the  location 
x  =  0,  y  —  r.  It  is  interesting  to  consider  the  behaviour  of  total  analysis  times  as  a 
function  of  the  number  of  degrees  of  freedom  (DOFs).  This  behaviour  is  plotted  in 
Figure  3  using  a  log-log  scale.  While  the  responses  are  different,  it  is  seen  that  analysis 
times  are  roughly  the  same  for  p-elements  and  h-elements  for  the  same  numbers  of 
degrees  of  freedom.  The  differences  in  analysis  times  are  further  discussed  in  Section 

4.2  for  3D  analyses. 


Table  1:  PAFEC  results  for  2D  h-element  analyses. 


No. 

elements 

No. 

DOFs 

Total  analysis 
time  (seconds) 

Maximum  stress 

concentration  factor 

Error  (Ski) 

10 

77 

1 

2.67700 

-16.80% 

40 

275 

2 

3.19310 

-0.76% 

160 

1031 

6 

3.26483 

1.47% 

640 

3983 

24 

3.24225 

0.77% 

2560 

15647 

142 

3.22587 

0.26% 

10240 

62015 

1173 

3.21923 

0.05% 

2  MSC.NASTRAN  outputs  an  equivalent  number  of  active  DOFs  for  p-element  analyses. 
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Table  2:  MSC.NASTRAN  results  for  2D  p-element  analyses. 


P-element 

order 

No.  Total  analysis 
DOFs  time  (seconds) 

Maximum  stress 

concentration 

factor 

Error  (Ski) 

MSC.NASTRAN 

error  estimate 

1 

588 

5.6 

2.93351 

-8.83% 

6.10% 

2 

744 

5.6 

3.28906 

2.22% 

10.05% 

3 

1356 

6.6 

3.23722 

0.61% 

1.05% 

4 

2608 

9.6 

3.22598 

0.26% 

0.62% 

5 

3660 

14.6 

3.21738 

-0.006% 

0.09% 

6 

4912 

25.4 

3.21742 

-0.004% 

0.04% 

7 

6364 

40.2 

3.21756(Kt(exact)) 

0% 

0.03% 

It  can  be  seen  that  there  was  very  good  agreement  in  the  final  converged  values  of 
stress  concentration  factor  for  both  element  types.  Here  the  h-element  analyses 
converged  to  a  maximum  stress  concentration  factor  of  3.21923  with  the  difference 
between  the  last  two  analyses  being  0.206%.  Here,  as  indicated  in  Table  1,  the  mesh 
density  was  increased  until  the  convergence  behaviour  was  monotonic  for  3  successive 
analyses.  The  p-element  analyses  converged  to  a  maximum  stress  concentration  factor 
of  3.21756,  with  the  difference  between  the  last  two  analyses  (orders  six  and  seven) 
being  only  0.004 %3.  Hence  for  the  p-element  analyses,  the  convergence  behaviour  can 
be  considered  as  being  very  good,  with  only  small  oscillation  about  the  final  converged 
value  as  the  element  order  was  increased.  This  is  important  to  note  since  it  indicates 
good  robustness  in  the  approach. 

Figure  4  shows  the  relationship  between  the  absolute  value  of  the  percentage  error  and 
the  number  degrees  of  freedom  for  both  h-element  and  p-element  cases,  plotted  on  a 
log-log  scale.  The  lines  of  best  fit  indicate  that  the  rate  of  convergence  (or  decrease  in 
error)  is  significantly  faster  for  p-elements.  The  results  in  Table  2  also  show  that  the 
MSC.NASTRAN  error  estimator  gives  a  useful  and  conservative  estimate  of  the  'true' 
error,  for  all  2D  p-element  orders  greater  than  one. 

3.3  Empirical  2D  'P-Element  Order'  and  'H-Element  Mesh  Density' 
Relationship 

Here  a  relationship  is  sought  between  the  order  of  a  2D  p-element  model  and  the  mesh 
density  of  an  h-element  model  with  the  same  mesh  discretization  error,  based  on  the 
results  given  in  the  preceding  section,  for  an  approximately  uniform  mesh  density. 


3  Based  on  the  discussion  in  Section  2,  the  'exact1  solution  for  the  2D  analyses  is  taken  as  Kt(exact) 
=  3.21756. 
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This  is  sought  on  the  basis  that  the  relationship  may  be  applicable  to  specifying  a  high 
quality  non-uniform  h-element  mesh  (in  one  iteration)  based  on  the  results  of  an 
automated  p-element  convergence  study.  For  example  an  accurate,  high  order 
p-element  analysis  using  a  relatively  course  geometric  mesh  in  MSC.NASTRAN  may 
be  used  to  guide  the  design  of  a  high  quality  h-element  mesh  in  ABAQUS  (which  does 
not  support  p-elements),  without  the  need  for  time  consuming  multiple  geometric 
mesh  refinement  iterations. 

First  we  can  postulate  that  the  negative  log  of  the  error  is  a  linear  function  of  p-element 
order.  Figure  5  represents  the  data  in  this  form,  where  the  fitted  equation  is  also  given4. 
The  correlation  coefficient  (R2)  is  0.9453,  which  indicates  that  the  fit  of  the  data  to  the 
hypothesised  relationship  is  very  good.  Next  for  h-elements,  it  is  hypothesised  that  the 
negative  log  of  the  error  has  a  power  dependence  on  the  number  of  h-element 
divisions.  Figure  6  shows  the  line  of  best  power  fit  and  the  equation  of  that  line.  In  this 
case,  the  correlation  coefficient  (R2)  is  0.7419.  The  poorer  correlation  is  due  to  the 
greater  oscillation  of  results  for  low  numbers  of  h-element  divisions.  We  can  now 
equate  the  response  equations  given  in  Figure  5  and  Figure  6  to  obtain: 

h  =  (0.51 12/7  -I-  0.1586)3  0469  (2) 

where  p  is  the  order  of  the  p-element  and  h  is  the  number  of  h-element  divisions  per 
p-element  division.  Equation  (2)  is  represented  graphically  in  Figure  7. 

It  is  important  to  note  that  this  equation  is  based  on  FE  results  for  models  with 
approximately  uniform  h-element  mesh  densities,  and/  or  uniform  p-element  order 
distributions.  Potential  applicability  to  non-uniform  cases  is  discussed  subsequently. 


3.4  Trial  of  H-P  Element  Density  Relationship  for  Non-Uniform 
Meshes 

In  this  section  the  empirical  relationship  obtained  in  the  previous  section,  based  on 
approximately  uniformly  distributed  mesh  density  (i.e.  h-element  size  or  p-element 
order),  is  assessed  for  the  case  where  the  mesh  density  is  significantly  non-uniform.  For 
this  trial  the  base  p-element  FE  model  had  40  elements,  and  is  shown  in  Figure  2(b). 
The  steps  used  to  conduct  the  trial  were  as  follows;  (i)  analyse  die  base  p-element 
model  with  the  auto-order  option,  (ii)  apply  the  h-p  relationship  from  Equation  (2)  to 
create  an  equivalent  h-element  model  using  the  auto-order  results,  and  (iii)  analyse  the 
h-element  model,  and  compare  the  error  from  the  p-element  and  h-element  models. 

The  auto-order  p-element  analysis  starting  p-order  was  set  at  one,  and  the  maximum 
allowed  p-order  was  five.  The  iterative  analysis  was  terminated  when  the  maximum 
node  stress  discontinuity,  from  adjacent  elements,  became  less  than  one  per  cent. 


4  Microsoft  Excel  97  SR-1  was  used  to  determine  the  line  of  best  fit. 
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Figure  8  shows  the  final  p-element  orders  for  each  element,  from  the  MSC.NASTRAN 
analysis,  which  range  from  two  to  five.  Table  3  shows  the  number  of  edge  divisions 
(rounded  to  the  nearest  integer)  used  in  the  equivalent  h-element  model  as  a  function 
of  the  p-element  order.  If  elements  sharing  a  common  edge  had  different  orders,  the 
higher  order  was  used  to  divide  the  edge.  Figure  9  shows  the  resulting  h-element 
model,  which  contained  6492  elements.  It  is  noteworthy  that  Figure  9  indicates  that  the 
maximum  density  of  h-elements  radiates  from  the  45-degree  positions  where  the  shear 
stresses  are  the  greatest.  This  confirms  that  typically  for  an  approximately  uniform 
distribution  of  h-element  sizes,  FE  errors  are  greatest  where  element  shear  stresses  are 
highest.  This  feature  appears  not  to  be  well  known,  as  some  analysts  incorrectly 
assume  that  the  greatest  mesh  density  should  be  at  the  location  of  peak  direct  stress. 

Table  3:  H-element  edge  divisions  (to  nearest  integer)  equivalent  to  p-element  order. 


P-element 

Equivalent  h-element 

order 

edge  divisions 

2 

2 

3 

5 

4 

11 

5 

21 

Results  of  the  two  analyses  are  compared  using  the  error  measure  defined  in  Section  2, 
where  Kt(exact)  is  given  in  Section  3.2.  The  p-element  analysis  gave  an  error  in  Kt  of 
0.555%,  whereas  the  h-element  analysis  error  was  0.087%.  Hence  the  empirical 
relationship  has  been  very  useful,  providing  a  conservative  h-element  analysis  result, 
which  has  less  error  than  the  initial  accurate  p-element  analysis.  If  it  is  desired  to 
reduce  the  level  of  conservatism  it  may  be  useful  to  use  the  average  of  the  two  element 
orders  to  divide  the  shared  edge,  as  opposed  to  using  the  higher  p-element  order. 

Hence  in  summary,  the  results  obtained  in  this  section  indicate  that  assuming  a  2D 
p-element  analysis  with  automatic  ordering  has  been  performed,  the  resultant  element 
orders,  in  conjunction  with  equation  (2),  can  guide  the  modeller  in  efficiently  designing 
a  h-element  mesh  for  further  complex  analysis  for  which  p-element  technology  is  not 
available. 
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4.  Convergence  of  3D  P-Element  and  H-Element 

Benchmarks 


4.1  Convergence  Analyses 

This  section  reports  on  analyses  with  several  variations  of  3D  p-element  and  h-element 
meshes.  The  following  specific  element  cases  were  analysed,  (i)  MSC.NASTRAN 
8-noded  brick  p-elements,  (ii)  MSC.NASTRAN  8-noded  brick  h-elements, 
(iii)  MSC.NASTRAN  20-noded  iso-parametric  brick  h-elements,  and  (iv)  PAFEC 
20-noded  iso-parametric  brick  h-elements.  A  model  composed  of  20  elements,  as 
shown  in  Figure  10  (a),  was  the  starting  point  for  all  sets  of  h-element  analyses.  In  the 
case  of  the  MSC.NASTRAN  8-noded  and  20-noded  brick  h-elements,  the  mesh  density 
was  uniformly  varied  in  three  dimensions  from  one  element  through  the  half  thickness 
of  the  plate  to  eight  elements  through  the  half  thickness.  Typical  meshes  are  shown  in 
Figure  10  (a-c).  Similarly,  for  the  PAFEC  20-noded  brick  h-elements,  the  mesh  density 
was  varied  from  one  to  six  elements  through  the  plate  half  thickness.  The  model  for  the 
p-element  analyses  was  the  same  as  the  base  model  for  the  h-element  analyses,  with 
the  polynomial  order  of  the  elements  varying  from  one  through  to  eight. 

4.2  Convergence  Results 

The  key  output  results  obtained  for  the  various  analysis  cases  are  given  in  Table  4  to 
Table  7,  as  a  function  of  mesh  size  or  element  order.  Here  to  determine  the  error,  the 
exact  value  was  taken  as  Kt(exact)  =  3.3754  (from  Table  10,  Section  5.2). 


Table  4:  MSC.NASTRAN  results  for  8-noded  brick  p-element  analyses. 


P-element 

order 

No. 

DOFs 

Total  analysis  Maximum  stress 
time  (seconds)  concentration  factor 

Error  (SKt) 

1 

351 

4.8 

3.1494 

-6.70% 

2 

462 

5.3 

3.4975 

3.62% 

3 

1053 

7.0 

3.4413 

1.95% 

4 

2684 

12.9 

3.3983 

0.68% 

5 

3895 

30.0 

3.3867 

0.34% 

6 

5346 

50.6 

3.3739 

-0.04% 

7 

7217 

114.8 

3.3775 

0.06% 

8 

9568 

267.8 

3.3762 

0.02% 
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Table  5:  MSC.NASTRAN  results  for  8-noded  brick  h-element  analyses. 


No. 

elements 

No. 

DOFs 

Total  analysis 
time  (seconds) 

Maximum  stress 

concentration  factor 

Error(SKt) 

20 

131 

3.4 

2.6863 

-20.42% 

160 

742 

5.4 

3.1533 

-6.58% 

540 

2193 

11.1 

3.3084 

-1.98% 

1280 

4844 

21.7 

3.3651 

-0.30% 

2500 

9055 

41.6 

3.3901 

0.43% 

4320 

15186 

77.2 

3.4015 

0.77% 

6860 

23597 

139.9 

3.4072 

0.94% 

10240 

34648 

224.3 

3.4097 

1.02% 

Table  6:  MSC.NASTRAN  results  for  20-noded  brick  h-element  analyses. 


No. 

elements 

No. 

DOFs 

Total  analysis 
time  (seconds) 

Maximum  stress 

concentration  factor 

Error(Sict) 

20 

442 

4.2 

3.4646 

2.64% 

160 

2684 

12.7 

3.5352 

4.74% 

540 

8166 

38.5 

3.4792 

3.07% 

1280 

18328 

98.6 

3.4463 

2.10% 

2500 

34610 

226.6 

3.4265 

1.52% 

4320 

58452 

627.8 

3.4144 

1.15% 

6860 

91294 

N.A5 

3.4044 

0.86% 

10240 

134576 

N.A 

3.3980 

0.67% 

5  The  final  two  analyses  were  run  on  a  different  hardware  platform,  and  hence  the  total  analysis 
times  cannot  be  validly  compared  with  other  results. 
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Table  7:  PAFEC  results  for  20-noded  brick  h-element  analyses. 


No. 

elements 

No. 

DOFs 

Total  analysis 
time  (seconds) 

Maximum  stress 

concentration  factor 

Error(SKt) 

20 

442 

3.8 

3.4187 

1.28% 

160 

2684 

30.3 

3.5212 

4.32% 

540 

8166 

173.2 

3.4603 

2.52% 

1280 

18328 

992.5 

3.4330 

1.71% 

2500 

34610 

4052.8 

3.4166 

1.22% 

4320 

58452 

13239.4 

3.4062 

0.91% 

The  results  in  Table  5  and  Table  6  show  that  the  MSC.NASTRAN  h-element  analyses 
with  8-noded  bricks  required  about  the  same  analysis  time  per  degree  of  freedom  as 
compared  to  the  MSC.NASTRAN  h-element  analyses  with  20-noded  bricks. 
Comparison  of  Table  6  and  Table  7  shows  that  unlike  the  2D  cases,  the  3D  PAFEC 
analysis  times  for  h-elements  was  about  an  order  of  magnitude  greater  than  for  the 
MSC.NASTRAN  h-element  analysis  times  for  the  larger  models.  Although  the  PAFEC 
time  includes  time  for  mesh  generation  (unlike  MSC.NASTRAN),  it  is  believed  that 
most  of  the  increase  is  due  to  other  differences  between  the  two  codes  and/or  how  the 
specific  features  of  the  codes  interact  with  computer  system  configuration.  Table  4 
shows  that  the  times  per  degree  of  freedom  for  the  p-element  analyses  were  similar  to 
the  MSC.NASTRAN  h-element  runs  for  lower  element  orders  (one  to  four),  but  were 
greater  for  higher  element  orders  (five  to  eight).  However,  one  key  point  from  Tables  4 
and  6,  is  that  the  total  solution  time  for  a  uniform  order  p-element  analysis  is  much  less 
than  that  for  uniform  density  h-element  anaysis,  for  a  typical  error  of  1  %. 

For  convenience  the  maximum  stress  concentration  factor  results  in  Table  4  to  Table  7 
are  summarised  in  Figure  11.  It  is  interesting  to  note  that  for  all  element  types  there  is 
good  agreement  in  Kt  to  within  1%  error.  Here  it  appears  all  cases  are  converging 
monotonically  to  the  true  Kt  result,  except  the  MSC.NASTRAN  8-noded  brick 
h-elements,  which  are  converging  to  a  result  approximately  1%  higher  than  the  true 
value.  As  expected  the  p-element  analyses  showed  the  highest  degree  of  convergence, 
reaching  to  within  0.06  percent  and  0.02  percent  of  the  exact  solution  with  p-elements 
of  order  seven  and  eight  respectively,  and  with  the  least  degrees  of  freedom. 

Figure  12  gives  the  response  of  the  absolute  value  of  the  percentage  error  as  a  function 
of  the  number  of  degrees  of  freedom,  for  the  four  sets  of  analyses.  As  expected  it  can  be 
seen  that  the  p-element  method  is  efficient  and  accurate  since  it  requires  significantly 
less  analysis  time  to  achieve  a  solution  within  one  percent  of  the  'exact'  solution.  By 
contrast,  both  20-noded  brick  methods  require  large  numbers  of  degrees  of  freedom 
and  analysis  times. 
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4.3  Empirical  3D  'P-Element  Order'  and  'H-Element  Mesh  Density' 
Relationship 

In  this  section  an  empirical  relationship  is  derived  between  the  order  of  a  3D  p-element 
model  with  uniform  order  distribution,  and  the  uniform  mesh  density  of  an  h-element 
model  (20-noded  brick)  with  the  same  mesh  discretization  error.  Here  a  similar  method 
was  used  as  described  for  the  2D  case  in  Section  3.3.  Initially  it  is  hypothesised  that  the 
negative  log  of  the  error  is  a  linear  function  of  p-element  order,  where  Figure  13 
represents  the  data  in  this  form.  Here  the  line  of  best  linear  fit6  is  given,  which  has  a 
correlation  coefficient  (R2)  of  0.9635.  For  the  20-noded  h-elements,  it  is  hypothesised 
that  the  negative  log  of  the  error  has  a  linear  dependence  on  the  log  of  the  number  of 
h-element  divisions.  Using  this  assumption,  Figure  14  shows  the  lines  of  best  fit  for 
both  the  MSC.NASTRAN  and  PAFEC  cases,  where  the  correlation  coefficients  (R2)  are 
0.9891  and  0.9986  respectively7.  The  three  correlation  coefficients  all  indicate  a  very 
good  fit  of  the  data  to  the  hypothesised  relationships.  Hence  combining  the  fitting 
equations  given  in  Figure  13  and  Figure  14  results  in: 

^  _  ^(0.6076^-0.2267) 

for  MSC.NASTRAN  h-elements,  and 

£  _  e(0.6096/,-0.349l)  ^ 

for  PAFEC  h-elements,  where  p  is  the  order  of  the  p-element  and  h  is  the  number  of 
h-element  divisions  per  p-element  division.  Equation  (3)  and  Equation  (4)  are 
represented  graphically  in  Figure  15.  It  is  noted  that  comparison  of  Figure  7  and  Figure 
15  shows  that  the  h-p  element  2D  and  3D  relationships  yield  about  the  same  values  for 
h-element  mesh  density  for  a  given  p-element  order. 

It  is  important  to  note  that  the  above  empirical  3D  relationship  has  been  trailed  in  a 
practical  application  involving  notch  stress  evaluation  in  a  large-scale  FE  model  of  the 
F-lll  wing  pivot  fitting.  In  that  work  the  local  h-element  mesh,  around  critical  fuel 
flow  vent  holes,  was  designed  so  that  a  non-linear  analysis  could  be  performed  using 
ABAQUS,  by  undertaking  an  initial  p-element  analysis  of  a  relatively  coarse  mesh  in 
MSC.NASTRAN.  This  proved  to  be  useful  in  providing  accurate  complex  h-element 
mesh  with  minimal  iterations  [13].  Here  it  should  be  noted  that  p-elements  are 
currently  not  available  for  non-linear  use,  in  any  commercial  code. 


6  Microsoft  Excel  97  SR-1  was  used  to  determine  the  line  of  best  fit. 

7  Results  for  both  MSC.NASTRAN  and  PAFEC  analyses  with  one  h-element  division  do  not 
follow  the  trend,  and  hence  are  excluded  from  curve  fitting. 
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5.  P-Element  Modelling  of  Through  Thickness 

Stresses 


5.1  Analysis 

In  this  section  further  work  is  reported  which  focuses  solely  on  several  variations  of  3D 
p-element  analyses.  Here  the  aim  is  to  investigate  in  detail  the  accuracy  of  different 
combinations  of  mesh  densities  and  element  orders  for  determining  the  peak  stress 
distribution  through  tire  thickness  of  a  plate,  at  the  hole  boundary.  In  particular,  to 
determine  whether  the  known  non-uniform  through  thickness  variation  in  Kt  can  be 
accurately  modelled  with  a  single  element.  As  indicated  in  Section  2,  here  32-noded 
brick  p-elements  were  used,  so  that  two  midside  nodes  could  be  used  to  define  the 
curvature  of  elements  around  the  hole  perimeter.  Table  8  summarises  the  four 
p-element  FE  geometric  models  used  in  this  investigation,  and  Figure  16  shows  the 
corresponding  FE  meshes.  It  should  be  noted  that  the  first  listed  mesh  was  a  full  plate 
thickness  model  while  the  other  three  were  half  thickness  models.  The  results  obtained 
are  to  be  subsequently  compared  to  results  from  an  accurate  PAFEC  analysis  with  a 
highly  refined  non-uniform  h-element  mesh  [10].  For  reference,  the  local  detail  of  the 
mesh  for  the  PAFEC  analysis  is  shown  in  Figure  17.  For  each  of  the  present  analyses, 
the  p-element  order  (uniform  throughout  the  mesh)  was  increased  by  one  until  the 
change  in  the  maximum  stress  concentration  factors  between  successive  analyses  was 
less  than  0.05%. 

Table  8:  P-element  FE  modelling  for  through  thickness  stress  analysis. 


Total  no.  of 
elements  in  model 

Equivalent  no.  of  elements 
through  full  thickness 

Symmetry 

used 

20 

1 

2  planes 

20 

2 

3  planes 

160 

4 

3  planes 

1280 

8 

3  planes 

5.2  Results  and  Discussion 

The  maximum  stress  concentration  factor  results  at  various  locations  through  the  plate 
thickness  are  given  for  each  of  the  four  meshes  as  a  function  of  p-element  order  in 
Figures  18  to  24.  These  are  sequenced  as  follows:  (i)  one  element  through  the  full  plate 
thickness  in  Figures  18  and  19,  (ii)  two  elements  through  the  full  plate  thickness  in 
Figures  20  and  21,  (iii)  four  elements  through  the  full  plate  thickness  in  Figures  22  and 
23,  and  (iv)  eight  elements  through  the  full  plate  thickness  in  Figure  24.  In  each  of  the 
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four  mesh  cases,  very  good  converged  stress  results  were  achieved  as  the  p-element 
order  was  increased.  As  can  be  expected,  more  pronounced  oscillation  about  the  exact 
stress  distribution  occurs  for  low  p-element  orders  for  the  least  refined  mesh,  as  shown 
in  Figure  18.  Conversely,  for  the  more  refined  meshes,  there  is  minimal  oscillation  even 
for  low  p-element  orders,  as  evidenced  in  Figure  24.  To  allow  more  precise 
comparison,  the  convergence  results  as  a  function  of  p-element  order  are  given  in  Table 
9,  where  the  difference  in  the  maximum  stress  concentration  factors  between 
successive  analyses  was  less  than  0.05%.  It  can  be  seen  that  the  four  p-element  models 
converged  to  maximum  stress  concentration  factors  in  the  range  3.3754  to  3.37568.  This 
compared  to  a  maximum  stress  concentration  factor  of  3.3790  from  the  PAFEC  analysis 
with  the  highly  refined  non-uniform  h-element  mesh  [10].  Similarly,  Table  10  shows 
the  p-element  order  required  for  each  model  to  obtain  accurate  (but  less  stringent) 
results  to  better  than  0.5%  error. 


Table  9:  P-element  order  required  for  converged  results  better  than  0.05%  difference  between 
successive  iterations  for  different  element  meshes. 


Equivalent  no.  of  elements 
through  full  thickness 

P-element 

order 

Converged  maximum  stress 
Concentration  factor 

1 

8 

3.3756 

2 

8 

3.3756 

4 

6 

3.3754 

8 

5 

3.3755 

Table  10:  P-element  orders  required  for  accurate  results  with  error  less  than  0.5%  for  different 
element  mesltes. 


Equivalent  no.  of  elements 
through  full  thickness 

P-element 

order 

Error(SKt) 

MSC.NASTRAN 

error  estimate 

1 

5 

0.30% 

0.014% 

2 

5 

-0.33% 

0.0033% 

4 

4 

-0.11% 

0.00061% 

8 

3 

-0.15% 

0.026% 

8  The  smallest  difference  between  consecutive  results  occurred  for  four  elements  through  the 
plate  thickness  with  p-element  orders  5  and  6.  Hence  3.3754  is  considered  to  be  the  'exact' 
solution  for  analyses  with  3D  (brick)  p-elements. 


13 


DSTO-TR-1151 


An  assessment  can  also  be  made  of  the  accuracy  of  the  MSC.NASTRAN  grid  point 
error  estimator9  for  3D  p-elements.  From  Table  10  it  can  be  seen  that  unlike  the  2D 
analysis  results,  the  3D  MSC.NASTRAN  error  estimator  underestimates  the  true  error 
for  these  higher  orders.  However,  it  is  worth  noting  that  in  almost  all  cases  with 
p-element  order  greater  than  two,  the  error  in  maximum  stress  concentration  factor  is 
less  than  one  percent. 

In  summary,  an  important  result  is  that  unlike  h-elements,  p-elements  (of  order  five 
and  above)  can  be  used  for  the  present  loading  to  predict  accurate  through-thickness 
peak  stresses  in  plates  modelled  with  only  one  element  through  the  full  plate  thickness. 
Based  on  the  present  results,  this  feature  has  been  used  in  some  recent  large-scale  FE  of 
RAAF  airframe  components,  to  save  significant  resources  in  FE  mesh  generation  time, 
and  numerical  computation  [13-17]. 


6.  Conclusion 


For  a  typical  local  stress  concentrating  structural  feature,  namely  a  circular  hole  in  a 
remotely  loaded  plate,  both  2D  and  3D  analyses  have  been  undertaken  to  assess  the 
accuracy  and  computational  resources  needed  for  various  uniformly  distributed 
h-element  and  p-element  mesh  densities,  and  p-element  orders.  Very  accurate 
converged  results  were  obtained  for  both  element  formulations.  However,  as  expected, 
p-elements  offered  large  reductions  in  analysis  times.  Using  these  results,  an  empirical 
relationship  was  developed  for  both  2D  and  3D  meshes  to  determine  the  equivalent 
non-uniform  h-element  mesh  density  which  will  yield  approximately  the  same 
accuracy  as  a  given  non-uniform  order  and  converged  p-element  mesh.  Through  a 
numerical  example,  the  veracity  of  this  approach  was  confirmed,  and  hence  this 
approach  provides  a  valuable  transferable  tool  in  designing  cost-effective  h-element 
meshes,  obviating  the  need  for  multiple  mesh  refinement  iterations.  This  is  particularly 
important  for  large-scale  finite  element  modelling  where  p-element  technology  is  not 
available  currently,  such  as  the  ABAQUS  and  PAFEC  codes,  which  are  used  in  AED  for 
non-linear  analysis  and  automated  shape  optimisation  respectively.  Finally,  an 
investigation  was  undertaken  which  demonstrated  that  p-elements  of  order  five  and 
above  can  be  used  to  predict  accurate  peak  through-thickness  stresses  in  plates 
modelled  with  only  one  element  through  the  wall  thickness. 


9  The  MSC.NASTRAN  error  estimator  is  described  in  Section  3.2. 
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(a)  One  division 
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(d)  Eight  divisions 


(e)  Sixteen  divisions 


Figure  2:  Two-dimensional  finite  element  meshes  for  h-element 
divisions  per  half  width  of  plate  specimen. 


(c)  Four  divisions 
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Figure  3:  Relationship  between  analysis  time  and  degrees  of  freedom  for  2D  h-  and  p-element 
models. 


No.  degrees  of  freedom 

Figure  4:  Relationship  between  absolute  value  of  percentage  error  and  degrees  of  freedom  for  2D 
h-  and  p-element  models. 
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Figure  5:  Relationship  betiveen  log(abs.  error)  and  p-element  order  for  2D  models. 
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Figure  6:  Relationship  between  log(abs.  error)  and  h-element  divisions  for  2D  models. 
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Equivalent  h-element  divisions 


(a)  Complete  model 


(b)  Local  detail  near  hole 


Figure  9:  Two-dimensional  h-element  FE  model  based  on  empirical  h-p  relationship  and  final 
p-element  orders. 


(a)  One  element  through  half  thickness 


(b)  Two  elements  through  half  thickness 


(c)  Four  elements  through  half  thickness 


Figure  10:  Typical  3D  finite  element  meshes  for  h-element  analyses  of  plate  specimen. 
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Figure  11:  Variation  of  maximum  stress  concentration  factor  Kt  with  degrees  of  freedom  for  3D 
h-element  and  p-element  models. 
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Figure  12:  Variation  of  absolute  percentage  error  with  degrees  of  freedom  for  3D  h-element  and 
p-element  models. 
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Figure  13:  Relationship  between  ln(abs.  error)  and  p-element  order  for  3D  element  models. 


Number  of  h-element  divisions  per  p-element  edge 


Figure  14:  Relationship  between  ln(abs.  error)  and  h-element  divisions  for  3D  element  models. 


DSTO-TR-1151 


€  g  40 


E  o  30 
0  £ 

0  _o 

A  0 _ 


P-element  order 

Figure  15:  Relationship  between  p-element  order  and  equivalent  h-element  divisions  for  3D  models. 


(a)  One  element  through  full  thickness 


(b)  One  element  through  half  thickness 


(c)  Two  elements  through  half  thickness  (d)  Four  elements  through  half  thickness 

Figure  16:  P-element  3D  FE  models  for  through  thickness  stress  analysis. 
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Figure  18.  Stress  concentration  factor  Kt  along  bore  of  hole  for  analysis  with  one  element 
through  full  plate  thickness  and  p-element  orders  1  to  4. 
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Figure  19:  Stress  concentration  factor  Kt  along  bore  of  hole  for  analysis  with  one  element 
through  full  plate  thickness  and  p-element  orders  5  and  6. 
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Figure  20:  Stress  concentration  factor  Kt  along  bore  of  hole  for  analysis  with  two  elements 
through  full  plate  thickness  and  p-element  orders  1  to  4. 


Figure  21:  Stress  concentration  factor  Kt  along  bore  of  hole  for  analysis  with  two  elements 
through  full  plate  thickness  and  p-element  orders  5  and  6. 
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Figure  22:  Stress  concentration  factor  Kt  along  bore  of  hole  for  analysis  with  four  elements 
through  full  plate  thickness  and  p-element  orders  1  to  4. 
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Figure  23:  Stress  concentration  factor  Kt  along  bore  of  hole  for  analysis  with  four  elements 
through  full  plate  thickness  and  p-element  orders  5  and  6. 
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Normalised  distance  along  bore  of  hole  (z/t) 

Figure  24:  Stress  concentration  factor  Kt  along  bore  of  hole  for  analysis  with  eight  elements 
through  fall  plate  thickness  and  p-element  orders  1  to  4. 
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